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A numerical simulation for "running wet" aircraft anti-icing systems is developed. The model 
includes breakup of the water film, which exists in regions of direct impingement, into individual 
rivulets. The wetness factor distribution resulting from the film breakup and the rivulet configuration 
on the surface are predicted in the numerical solution procedure. The solid wall is modeled as a multi- 
layer structure and the anti -icing system used is of the thermal type utilizing hot air and/or electrical 
heaung elements embedded within the layers. Details of the calculation procedure and the methods used 
are presented. 


Nomenclature 

Cp = specific heat 

F = wetness factor 

hi = heat transfer coefficient between the hot air and 
the inner surface of the wall 

h w = heat transfer coefficient between the outer 
surface of the wall and the runback water 
hoo heat transfer coefficient between the free stream 
and the outer surface of the wall 
k = thermal conductivity 

L y = latent heat of vaporization of water 

LWC = liquid Water Content 

M = molecular mass, number of grids across film 

m = runback water mass flow rate 

m = rate of mass transfer per unit area 

P = static pressure 

Pr = Prandtl number 

q = rate of heat transfer 

J9 

q = rate of heat transfer per unit area 

<f = rate of heat generation per unit volume 

R = rivulet radius 

r = recovery factor 

Sc =s Schmidt number 

V = flowfield velocity 
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w 

= velocity in a rivulet in the flow direction 

T 

= temperature 

a 

= thermal diffusivity 

p 

= rivulet contact angle with the solid surface 

9 

= equivalent rectangular film thickness of a 


rivulet 

X 

= ratio of rivulet width to wetness factor, or 


distance between two adjacent surface 


streamlines 

n 

= droplet collection efficiency 

n 

= dynamic viscosity of water 

S 

= area correction factor for heat loss from a 


rectangular film to the ambient 

p 

= density 

r 

= surface shear force (friction) 

<p oo 

= free stream relative humidity 

Subscripts 

a 

= anti-ice air 

at 

= anti-ice 

€ 

= property at edge of the boundary layer 

evap 

= evaporation from outer cowl surface 

f 

- liquid film 

imp 

= impingement on outer cowl surface 

m 

= solid wall composed of several layers 

r 

= rivulet 

vap, v 

= vapor, saturated vapor 

w 

= wall or runback water 

OQ 

= at free stream conditions 


1 


Superscripts 

l = layer number in the composite wall 
n = Az step level (grid number in the z -directi on) 

I. Introduction 

The problem of aircraft icing has been the focus of 
study of many researchers for a number years. The 
detrimental effects of ice accretion on critical surfaces can 
jeopardize flight safety as well as the overall aircraft 
performance. Consequently* accurate modeling and 
extensive study of the icing process are necessary. Two 
general methods of ice protection have been developed: De- 
icing methods for the intermittent removal of ice buildup 
by destroying the bond between the ice and the surface, and 
anti-icing methods for the prevention of ice formation on 
critical surfaces such as engine nacelles. 

The availability of high-speed digital computers has 
favored the use of numerical techniques and the 
development of computer codes to design and analyze ice 
protection systems. It is felt that the latter can minimize 
the cost associated with the required experimental testing 
by providing a tool that is at least capable of predicting 
preliminary results. 

Most studies related to aircraft icing have been 
committed to the prediction of ice shapes and the 
determination of their detrimental effects on aerodynamic 
performance of the aircraft components. At this time, 
research in running wet anti-icing systems is quite basic, 
and runback is treated in a primitive manner. The NASA 
Lewis Research Center has been a major contributor in 
conducting and sponsoring studies related to computer 
modeling of aircraft icing processes as well as 
experimental testing in its Icing Research Tunnel (IRT). 
As a result, LEWICE [1], an ice accretion prediction code, 
was developed for unprotected airfoil surfaces. The 
approach used in the modeling consists of performing 
mass and energy balances on the surface water. The 
wetness factor issue is ignored and the runback water is 
assumed to wet the entire surface at a particular location. 
Consequently, the amount of required heat to anti-ice the 
surface is under-estimated. 

Several investigators have produced different versions 
of the LEWICE code in order to improve it. To name a 
few, Cebeci et. al [2] modified the flowfield calculation 
module of the code to avoid the problem of multiple 
stagnation points. Yamaguchi et. al [3] proposed a multi- 
zone roughness model: a zone of uniform water film in the 
stagnation region, a rough zone of stationary water beads, 
and lastly, a zone where surface water run back as 
rivulets [4]. The runback water was recently modeled by 
AI-Khalil et. al [5,6,7] by incorporating a rivulet model. 
This paper is intended to present the numerical calculation 
procedures used including the most recent improvements 
of the latter model. 

II. Mathematical Model 

The runback model introduced earlier is based on a 
two-dimensional mathematical formulation. The surface 
water and the solid structure temper^ ires vary across their 


thicknesses and in the flow direction along a streamline on 
the surface. S panwise temperature dependence is assumed 
to be negligible. However, the latter is accounted for by 
performing energy balances on control volumes whose 
span wise widths extend between two adjacent streamlines 
on the aircraft surface. 

11.1 Runback Water 

IL1J Hydrodynamics 

The rate of water impingement on aircraft surfaces, 
due to the existence of undisturbed supercooled liquid water 
droplets in clouds, is relatively small. This and 
aerodynamic forces result in a very shallow water film 
flowing over the skin surface. Consequently, the surface 
water behavior is controlled not only by aerodynamic and 
body forces, but also by surface tension forces and surface 
roughness. 

In the direct impingement regions, i.e., in the 
neighborhood of the stagnation point, the water tends to 
wet the entire surface due to incoming droplets and due to 
water running back from upstream locations. However, at 
or downstream of the impingement limit, the liquid film 
could become unstable due to surface tension forces that 
cause the surface water to coalesce into individual streams, 
referred to as rivulets, separated by dry areas. 

A detailed study on the hydrodynamics and a stability 
analysis of surface water was presented in Ref. [5]. For 
completeness, some of the essential features are presented 
here without further disc ussion. T he film/rivulet flow in 
the streamwise (z) direction is caused by a shear force 
acting at the liquid-air interface. The latter force is 
obtained from the results of the skin friction factor 
computed from viscous aerodynamic calculations of the 
flowfield. 

A rectangular film model was chosen to 
mathematically represent the heat transfer process in a 
rivulet as shown in Fig 1. This model was found 
appropriate to the current problem for various reasons 
discussed in Ref. [6]. The criteria used for the new 
runback water configuration are as follows: 

• The wetness factor is preserved, i.e., the widths of a 
rectangular film is equal to that of its corresponding 
rivulet 

• The law of mass conservation requires equal mass flow 
rates in a rivulet and its equivalent rectangular film. 
This criterion enables one to compute the Him thickness 
5 f . 

• Mass loss due to evaporation is associated with a 
decrease in the rivulet size, i.e., its radius and, 
consequently, its base width that is also equal to the 
rectangular film width. This criterion enables one to 
update the value of the wetness factor. 

The velocity distributions within the film and the 
rivulet were derived and used to obtain the mass flow rates 
in each [5], as shown, respectively: 
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( 2 ) 

The second criterion is used, equating the above equations, 
to give: 


impinging droplets which are of value only in the direct 
impingement region. 

The rate of impingement per unit area, m"i mp , is 
calculated from the local value of the collection efficiency 
as shown: 



(3) 


where is a function of /3 derived in [7], The above 
equation shows the rectangular film thickness is directly 
proportional to its equivalent rivulet radius. This equation 
will later be used to update 5f when R is reduced due to 
evaporation. Procedures to determine the conditions and 
location for the breakup of the liquid layer flowing 
downstream of impingement regions were thoroughly 
discussed in Ref. [5]. The prediction of initial values of 
R and F at breakup was also described. 

11.1.2 Therm al Analysis: 


The principal objective of this study has been to 
analyze and predict the performance of anti-icing systems. 
In such applications, the worse case occurs at equilibrium 
state conditions. Consequently, the mathematical 
formulation of the heat transfer process is based on the 
steady-state energy equations. The unsteady equations are 
more relevant to de-icing applications. The runback water 
energy equation then follows: 


dz dy 2 


(4) 


where, y = J* w 

My) 

The above equation is based on the fact that 
conduction heat transfer within the liquid water in the z- 
direction (flow direction) is negligible compared to that in 
the y-direction (across film thickness). The solution of 
Eq. (4) requires two boundary conditions in the y- 
direction, one at the solid-liquid interface, and one at the 
liquid-air interface, and an initial condition (z=0). The 
latter condition requires knowledge of the water 
temperature at the stagnation point. Analytically, this is 
impossible because that depends on the final temperature 
solution in the water film and in the solid structure layers. 
However, this may be obtained numerically in an iterative 
procedure described in a later section. 

The boundary condition at the liquid-air interface is 
written as: 


k w ^= h-E 
dy 


T w -Te - 


2 C p< & j 


+ m, 


evap 


E, Lv 


- mimp Qv_ (To. - T w ) - (5) 

where the first and second term terms on the right-hand 
side represent heat loss to the ambient by convection and 
evaporation, respectively; and the third and fourth terms 
are the sensible and kinetic heat contributions of the 


m" imp = tj LWC V„ 


( 6 ) 


and the rate of evaporation per unit area is computed using 
the Chilton-Colbum heat-mass transfer analogy. This 
may be expressed as: 


fit tvap 


C P~ 


1 Pr j 273 Mh 7 o 

P v,w — P vap 

\ SC >air Moir 

P * ~ P v,w 


(7) 


where, 

P v , w = saturated vapor pressure at the local runback water 
temperature T w . 

P vap = local vapor pressure at the edge of the boundary 
layer at the local relative humidity. 

Application of Dalton's Law of partial pressures and 
knowledge of ambient conditions yields: 



where the relative humidity 0oo in a cloud is generally 
taken to be 100%. The saturation vapor pressure of water 
is written as function of temperature: 


Pv(T) = 2337 exp / 6789 T — 1 1 

l 1.293.15 T J 


-5.031 In 


293.15 


-) (9) 


where the units of P y and T are (Pa) ant (K), respectively. 

The recovery factor r in Eq. (5) accounts for viscous 
dissipation in the boundary layer and is approximated 
by [8]: 


r-l-Ji(l -PrZ,), 

vl 


f — laminar flow 

I 2 

— turbulent flow 
1 3 


( 10 ) 


The properties at the edge of the boundary layer, i.e., P e , 
T e , and V e are computed using the perfect gas relations for 
isentropic flow and the local values of the pressure 
coefficient obtained from a flowfield solver. 

Note that £, in Eq. (5), is an area correction factor to 
account for the area differences in the rivulet and the 
rectangular film models through which heat exchange with 
the ambient occurs. This factor is defined as the ratio of 
the rivulet free surface area to the upper surface area of the 
corresponding rectangular film. From geometric 
considerations, the following may be written: 

£ = JL (o<f < i) 

sinp 

S=i (r=i)j 
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This factor is less than 10% for contact angles smaller 
than 42°, and is unity for uniform film flow. 

The boundary conditions concerning the solid-liquid 
interface which represents the heat exchange between the 
solid wall and the runback water remains to be discussed. 
These two conditions will be presented with the energy 
equation of the wall structure since they are common 
between the two regions. 

11.2 Anti-ice Bleed Air: 

A widely used method of preventing ice formation is 
the hot air type due to its high reliability. In these 
systems, hot air is drawn from an intermediate or high 
stage compressor bleed port and ducted through passages to 
the leading edges of wings, empennages, engine nacelles 
or other critical areas. Due to the complexity of the flow 
of the anti-icing air inside the irregular duct shapes and the 
uniqueness of each design, a generalized model requires the 
following assumptions: 

1 . The heating requirement by such a system is generally 
specified by the amount of hot air supply, m a , and its 
delivery temperature at the stagnation region. 

2. The internal heat transfer coefficients hi, between the 
air and the inside surface of the structure, is assumed to 
be known from previous experience or from 
experimental testing on the particular system in 
consideration. 

3. The hot air temperature varies in the flow direction and 
is assumed to be lumped in the transverse direction. 

With regards to the above assumptions, the energy 
equation of the anti-ice air may be wriuen as: 

m a C Pa — = hi(z) [ T a (z ) - T m (y=0,z)] m) 


-k m ^ = q" a + h,(T a -T m ) 

dy 

where q" a i is an optional prescribed heat flux distribution. 
This value and h t may be set to zero for a perfectly 
insulated inner surface. The above equations were 
formulated as such to give the flexibility of modeling 
different systems. 

The two conditions that must be satisfied at each 
solid-solid interface between the wail layers are continuity 
of temperature and heat flux normal to the interface. As to 
the boundary conditions on the left side (stagnation point) 
and the right side of the wall, they may be extrapolated 
frorr^the solutions using insulated conditions, or they may 
be specified if the temperature distribution is known at 
either end. ~ 

The last boundary conditi ons jh at remain are those 
pertaining to the butermost layer at the solid-liquid 
(partially or fully wetted surface) and solid-air (dry surface) 
interfaces. They may be written as follows: 

T m =T w (0<F<1) (15) 

and, 

q^-k^F^- +h „ (1-F ) (7V7V- ±XL) (0£F< 1) (i 6 ) 
dy 2 Cpu, 

The first condition is only necessary in the fully or 
partially wetted regions. The second condition simply 
states that heat, q " m » is transferred from the wall 
proportionally through the wetted (to the water) and the 
dry (to the ambient) surface areas as defined by the wetness 
factor F. Note that T m , in Eq. (16), may be replaced with 
T w according to Eq. (15). 


Obviously, T a (z) depends on the solid wall temperature 
distribution which also depends on the runback water 
solutions. Therefore, the energy equations of those three 
regions must be solved simultaneously. 

11.3 Wall Structure: 

Based on the assumption that the wall temperature is 
dependent on the y and z-direction [6], the following 
energy equation may be written for each layer in the 
composite structure: 

A( z )^kl + ^Zk+il = 0 (13) 

. dz J dy 2 km 

where X(z) may be taken as the distance between two 
adjacent surface streamlines which make up the strip being 
analyzed. This distance is constant for a 2-dimensional 
flow over a surface. The above formulation allows one to 
model a heating element as one of the layers. If anti-icing 
is achieved by means of a hot air system alone, the value 
of q° may be conveniently set to zero for all the layers. 

The boundary condition at the inner surface of the wall 
may be written for the innermost layer as: 



III. Numerical Solution Techniques 
111.1 Runback Water: 


A fully implicit method was used to numerically 
solve Eq. (5) because of the positive stability 
properties [9]. Backwards differencing in the z-direction, 
and central differencing in the y-direction were employed. 
Applying this scheme to Eq. (5) and rearranging terms 
yields: 
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(17) 


for j=2,3,.. .^-1, where M is the total number of grid 
points across the film thickness (in the y-direction), and n 
is the grid number in the flow direction (z). Equation (17) 
is written for each corresponding node which results in a 
set of linear equations. The latter may be rewritten in a 
matrix form and solved using the Thomas Algorithm for 
tridiagonal system of equations. 

However, before carrying the solution, two equations, 
corresponding to 7=1 (solid-liquid boundary) and j=M 
(liquid-air boundary), are still required. A one-sided 
difference representation of Eqs. (16) and (5) is used for 
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this purpose, respectively. Equation (15) could have been 
used in the runback solution while Eq. (16) is used in the 
wall solution, instead. However, that procedure was found 
to be highly unstable. Thus: 


1 , h~(l -F)Ay 
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T J* I = 
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(18) 


where j= 1, and q" m is the rate of heat flux normal to the 
solid-liquid interface computed from the temperature 
distribution in the wall from the following: 


<?m — — km 


dTm 

dy 


( at solid-liquid interface ) 


(19) 


performed on a control volume of length Az t (distance 
between node n=l and n=2), thickness 5f, and a unit depth. 
Using Eq. (1) with A=1 and F=l, this yields: 

i^imp ~ ffltvap) A Z ] — ^ — Sf 

2/i 

Solving for fy gives: 

( ini “al) = V 2 pr^' (22) 

where t is taken as the average wall shear force between 
nodes n=l and n=2. 

The conservation of mass equation of the runback 
water may be readily obtained and shown to be: 


A second order finite differencing was used to compute 
the right-hand side of Eq. (19). At j=M y one may write, 
after rearranging terms: 


+ 


1 + (/,„$♦ 

k w 


m,„ p C Pw 


.) 


7 ' 1 

j - 


k w 
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( 20 ) 


m * +l = m* + X. n+1 Az (m-mp-m.'vap £ ) (23) 


Knowing the mass flow rate, the film thickness in the 
fully or partially welled regions may be derived from 
Eqs. (1). (2) and (3): 



2 p m"* 1 

p t X 


(^=1) (24) 


The above equations may now be solved for the 
temperatures at nodes ;'= 1 through M at location z+Az 
(i.e., n+1), knowing the nodal temperatures at location z 
(i.e., n). The evaporation term m" evap is computed using 
the temperature at z in order to preserve the linearity of the 
system of equations. 

The procedure described above requires knowledge of 
the water temperature at the stagnation point (z=0, or 
n=l). This is obtained by extrapolation from the 
temperature distributions at n=2 and n=3. Since the 
solution procedure is iterative, as discussed later, an initial 
guess is required to start the computations. This is 
achieved by performing mass and energy balances on a 
differential control volume of the surface water at the 
stagnation point, which yields the following approximate 
expression: 


rr 1 (initial) = 


+ h~ (Tc + 




imp — — — ) — ftlfyapFy + 


rVt 

2 C, 


Pcir J 


/< 


M imp Cp wm + ku» ) 


<7*1 


( 21 ) 


where q" m is estimated assuming that heat conduction 
within the solid structure occurs in the outward direction 
(y). Equation (21) is only used at the first iteration. In 
subsequent iterations, the extrapolation technique 
mentioned previously may be used. However, this caused 
slight fluctuations in the temperature distribution at the 
first few nodes (n= 1 ,2,3,4). The problem was remedied by 
setting the initial water temperature equal to the average 
temperatures of nodes n=2 and 3 without affecting the 
remaining results. 

In addition to an initial temperature, an initial water 
film thickness is required. A mass balance may be 


sr' 




P T Fi(P) 


fv 


^I(P) 
sin (3) 


(0<F<1) 


(25) 


where m r is the mass flow rate per rivulet. In the case 
where the runback water is flowing as rivulets (F<1), the 
wetness factor must be updated at each z-location. From 
geometric considerations, this is derived from 
F=(2f?sin^)/X where R is obtained from Eq. (3). 
However, if surface streamlines are not parallel (3-D flow), 
great care must be taken when evaluating X to account for 
variations in the distance between two surface streamlines 
which identify the strip being analyzed. 

The numerical solution of Eqs. (17), (18) and (20) 
requires the discretization of the water domain into grid 
points. Across the liquid layer thickness, equal spacing 
between the grid points was used. Along the flow 
direction, two zones were selected: direct impingement 
region, and downstream region. The grid spacing is 
constant in each zone, but is much smaller in the direct 
impingement region to accommodate for the rapidly 
changing variables due to the impinging water droplets and 
the flowfield characteristics. 

The current model was specifically developed for anti- 
icing applications where at least the minimum heat 
required to keep the surface water from freezing is supplied 
to the surface. This is because a two-dimensional phase- 
change model was found to be inappropriate since freezing 
will normally start at the liquid-air interface, which creates 
a problem in modeling the flow characteristics of the 
unfrozen water. However, since the temperature drop 
across the film thickness is small, the temperature may be 
assumed to be uniform across the layer. Therefore, when a 
freezing temperature, or lower, is obtained during the 
calculation process, an alternate method is used. This 
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consists of performing a macroscopic energy balance on 
the surface water to obtain the freezing fraction, such as 
done in the LEWICE code [1]. Nevertheless, the rivulet 
configuration and its prediction remain the same. This 
enables one to predict the amount and location of ice 
accumulation during a specified period of exposure time, 

111.2 Anti-ice Bleed Air: 


10. Outer totally/partially wetted surface node. 

11. Outer totally dry surface node. 

Energy balance equations for all node types are derived 
and presented below. The following definitions were used: 

/}=Ai = = and Q , = 

Ay ki Ayi Ayi (26) 


The governing energy equation of the anti-ice bleed 
air, Eq. (12), is a first order ordinary differential equation 
(ODE). Due to the arbitrary distribution of the heat 
transfer coefficient and the wall temperature at the inside 
surface of the solid structure, a numerical technique must 
be used to solve the latter equation. 

A forward finite difference scheme is only first order 
accurate, A more accurate and widely used technique for 
solving ODE's, is the fourth order Runge-Kutta 
method [10]. Knowing the temperature distribution in the 
wall, from the most recent iteration, the latter method is 
used to predict or update the hot air temperature 
distribution in the cowl. The result is subsequently used 
in the wall temperature solution at the next iteration. 

In cases where anti-icing is achieved by means other 
than the hot air type (i.e., m a = 0), the solution of Eq. (12) 
using the aforementioned technique should be avoided. 
Instead, the air in the cowl is considered to be stagnant and 
at a prescribed temperature. Also, when the internal heat 
transfer coefficient is zero (i.e., insulated inner surfaces), 
there is no need to solve Eq. (12) since the result is a 
constant air temperature which does not affect the wall 
temperature, and consequently the runback water 
temperature. 

III.3 Wall Structure: 


where / and /+1 indicate the layer numbers corresponding 
to a particular solid-solid interface. Note that in the 
following, node (ij) denotes the grid point at row and 
column n /\ and that A; represents the distance between the 
two surface streamlines, defining the width of the wall 
strip being analyzed, at column 

Node type 1; 
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= Azi\jL(hiT a +qai) + q° 
k [Ay 


Ti+\j-2 p TiV+i 
(28) 


A solution for the different layers in the wall structure 
may be obtained by direct approximation of the governing 
equation, Eq. (13), and the corresponding boundary 
conditions by finite differences. However, the control 
volume approach was chosen due to its accurate 
conservation properties [7]. Difference equations are 
derived by performing an energy balance on each control 
volume corresponding to a particular node. The control 
surfaces of each control volume are half way between the 
corresponding node and its adjacent surrounding nodes. 
There exist eleven types of nodes in the wall structure. 
These types are listed below and correspondingly numbered 
as shown in Fig 2. which illustrates a two-layer wall 
(note that the wall thickness dimension compared with its 
length is exaggerated for clarity): 

1. Totally internal node. 

2. Inner surface side node. 

3. Inner surface left-comer node, 

4. Inner surface right-comer node. 

5. Left-side internal node. 

6. Right-side internal node. 

7. Solid-solid interface internal node. 

8. Solid-solid interface left-side node. 

9. Solid-solid interface right-side node. 


Node type 3: 

If the temperature distribution is not specified at z=0, 
an insulated boundary condition is used: 
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Node type 4: 

Similarly, for unspecified temperature: 
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Node type 5: 

For unspecified temperature: 
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Node type 6: 

For unspecified temperature: 
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Node type 7: 
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Node type 8: 

For unspecified temperature: 
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Node type 9: 
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2A, 

~ + £ < i -rr-Tij+\=-Az-[q' > i + 
i2i 2 A; 

Node type 10: 

In this region, the node temperature is set equal to the 
local liquid temperature at the base of the film. This is 
achieved using a cubic spline interpolation technique 
because the interfacial grid points of the water and the 
solid do not coincide. 

Node type 1 1 : 

For consistency with the lower boundary of the liquid 
layer, a direct differencing of the equation representing 
convective heat loss to the ambient is applied. This 
yields: 

x . /\ h— dy . „ ho* Too 4y 

- TiM-i + (1 + J - ) T iM = L 

A A 


(36) 


III. 4 Solution Procedure: 

The required solutions are the temperature 
distributions in the anti-ice hot air, the solid structure, and 
the runback water. In addition, the surface water mass 
flow rate and the film/rivulets configuration must be 
determined. A simultaneous solution must be carried in 
the three regions: (1) runback water, (2) solid structure; 
and (3) anti-ice bleed air. This may not be accomplished 
in a single step due to the dependency of some boundary 
conditions of a particular region on the final solution in 
the adjacent region. This suggests the use of an iterative 
type of numerical solution between the three regions. 

The sequence of the steps utilized in the numerical 
solution iterative procedure may be listed as follows: 

(1) Estimate q m , in Eq. (18), at all nodes corresponding 

to the runback water at j= 1 . The procedure is to use 
a local one-dimensional heat transfer model from the 
wall to the free stream air (i.e., no conduction within 
the wall in the flow direction), assuming a fully dry 
surface. Any heat transfer generated due to electrical 
heating elements is assumed to flow outboard to the 
ambient. These assumptions were necessary to get 
the iterative solution started. 

(2) Compute the "initial" water temperature and the film 

thickness at the stagnation location from Eqs. (21) 
and (22), respectively. 

(3) Solve Eqs. (17), (18), and (20) for the runback water 
temperature distribution across the film thickness at 
the next z-location. Proceed with the solution of the 
latter equations by marching to the location of the 
impingement limit. This, of course, corresponds to 
the direct impingement region where the wetness 
factor is unity. Note that the runback water mass 
flow rate and the film thickness are updated using 
Eqs. (23) and (24), respectively, as the solution is 
brought to the next level. 

(4) From the impingement limit onward, check if the 
criteria for film breakup are met as the march 
proceeds downstream with the solution. If breakup 
occurs, the wetness factor and the rivulet 
configuration are predicted [5], Then proceed with 
the calculations for each step up to the end of the 
structure or up to the location where total 
evaporation occurs. The film thickness is updated 
using Eq. (24) or (25), and the wetness factor is 
updated by geometric considerations after each Az 
step. 

(5) Generally, a larger number of nodes is used in the 

runback water than in the wall at the solid-liquid 
interface. Thus, a cubic spline interpolation 
technique is used to predict the wall nodal 
temperatures, for node type 10, from the water nodal 
temperatures at the interface. 

(6) Setup the equations corresponding to convective 
boundary condition, Eq. (36), for nodes of type 11, 
if total evaporation occurs upstream of the end of the 
structure. 
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(7) Assume a constant anti-ice bleed air temperature equal 

to the delivery temperature at the stagnation region, 

(8) Setup the equations corresponding to the rest of the 

solid structure nodes, types 1 through 9. 

(9) Solve the linear system of equations for the wall 

nodal temperatures. This terminates the first 
iteration. 

(10) Compute q” m from Eq. (19) using the temperature 
distribution obtained in the previous step, and 
interpolate for the runback nodes using cubic splines. 
Under-relaxation of the latter values should be used 
to carry a stable solution as follows: 

n "n+\ _ n ”n a p ( n '*n + 1 _ n "n \ 

H m — H m ' ' \H m H m) 

where F is the under-relaxation factor and has a value 
between zero and unity. Its actual value depends on 
the particular problem under consideration, 

(11) Evidently, the solution would not converge in one 
iteration. Extrapolate for the initial water 
temperature as previously discussed, and compute the 
film thickness at stagnation from Eq. (22). 

(12) Repeat the runback water solution as described in 
steps (3) and (4). 

(13) Set the wall temperature at the solid-liquid interfacial 
nodes in the wet region by interpolation from the 
water solution of step (12). Also setup the 
convective boundary condition equations in the dry 
region as done in step (6). 

(14) Solve for the temperature of the anti-ice bleed air as 
described in section III. 2, using the most recent wall 
temperature distribution at the inner surface. 

(15) Setup the equations corresponding to the remaining 
wall nodes as in step (8), then solve the system of 
linear equations for the wall nodal temperatures. 

(16) Compare the solutions obtained in the previous step 
with the corresponding solutions of step (5). If the 
difference is within an acceptable tolerance, the 
solution is considered converged. Otherwise, 
perform another iteration by repeating the last few 
steps starting with step (9). 

IV. Sample Calculations and Discussion 

The primary purpose of this paper was to present the 
details of the mathematical development and the numerical 
solution techniques of the current model. Therefore, only 
one example problem will be considered in order to 
demonstrate the calculation procedure. However, several 
other cases were considered and presented in Ref. [11]. 
The complete solution to the problem is resolved in three 
major steps: (1) flowfield calculations, including the 
viscous layer near the wall; (2) individual water droplet 
trajectory calculations using the velocities calculated in the 
previous step; and, finally, (3) the heat transfer 
calculations for the anti-ice hot air, the solid structure, and 
the surface water. 

In the following example, the solid structure is 
assumed to be a NACA 0012 airfoil of chord length equal 
to 1.0 m, as illustrated in Fig. 3. The wall structure of 


the airfoil is composed of five layers, typical of some 
aircraft surfaces. Properties and dimensions of these layers 
are illustrated in Table 1. The electrical heater, center 
layer, is assumed to be turned off and heating of the 
surface is accomplished by spraying hot air on the inside 
surface of the cowling near the stagnation region. The air 
is delivered at a temperature of 200 °C and a mass flow 
rate of 0.1 Kg/sec per unit span wise distance. 

The ambient operating conditions are the following: 

• Flight Mach number = 0.25 

• Ambient static temperature = -12 °C 

• Ambient static pressure = 90.75 kPa 
■ Angle of attack = 0° 

• Cloud Liquid Water Content = 1.0 g/m 3 

• Relative humidity = 100% 

• Mean volume droplet diameter = 20 jim 

The flowfield around the airfoil was computed using 
the ARC2D code which solves the two-dim ensi onal thin 
layer Navier-Stokes equation. A hyperbolic grid generator 
was used to produce a C-type grid structure around the 
airfoil: 239 nodes along the surface and 55 nodes in the 
normal direction. Grids were packed near the wall for 
accurate prediction of the large gradients induced by 
viscous effects in these regions. The resulting pressure 
coefficient and friction coefficient distributions are 
illustrated in Figs. (4) and (5), respectively. These 
coefficients are defined as follows: 


C P = 


P -P- 

Ip- vl 

2 


and, C/= 


Ip. Vl 

2 


The first coefficient may be used to calculate properties at 
the edge of the boundary layer, and the second is used to 
compute the wall shear stress that cause the water to run 
back. 


A particle trajectory code was then used to produce a 
collection efficiency distribution on the surface, as 
illustrated in Fig. 6. Note that all the results presented 
thus far are symmetric between the upper and the lower 
surfaces of the airfoil. This is due to the fact that the flow 
angle of attack is zero and the airfoil geometry is 
symmetric. 

The final step involves the heat transfer calculations. 
The external convective heat transfer coefficients, between 
the wall surface and the ambient air, were computed using 
a sand roughness factor of ks/c=0.0002 [1]. The internal 
heat transfer coefficients, between the hot air and the inner 
surface of the afifoUcowl were arbitrarily assumed since 
they depend on the particular air jet nozzles design, the rate 
of air flow, and the air passages geometry. These 
coefficients are shown in Fig. 7. 

The procedure described earlier is applied, using the 
results thus far obtained, to solve for the problem variable 
parameters. The contact angle between the rivulets and the 
surface, when breakup of the film occurs, is assumed to be 
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/J=40°. The actual value of /? depends on the properties of 
the solid surface and its roughness. 

The resulting average temperature distribution of the 
anti-ice air inside the cowl is illustrated in Fig. 8. The air 
temperature drop across the entire length of the airfoil is 
approximately 85 °C. The drop occurs in a relatively 
smoother manner compared to that of the runback water 
average temperature, shown in Fig. 9. This is due to the 
distribution of the corresponding convective heat transfer 
coefficients. Since the solid wall conductivity is relatively 
larger than that of water, its average temperature 
distribution tends to be smoother as depicted in Fig. 9. 

The distribution of the heat flux leaving the outer 
surface of the airfoil is plotted in Fig. 10. The curve 
peaks are due to the peaks in the distribution of the 
external heat transfer coefficients which correspond to a 
transition from laminar to turbulent flow. Figures 1 1 and 
12 are plots of the runback water film thickness and the 
wetness factor, respectively. The sudden jumps in the 
curves correspond to the breakup of the uniform film in 
the direct impingement region (F=l) into individual 
rivulets (F<1). 

The corresponding distribution of the runback surface 
mass flow rate per unit spanwise distance is shown in 
Fig. 13. This system is clearly a running wet anti-icing 
system. Total evaporation may be better accomplished 
with electrical heating elements such that a large amount 
of heat is supplied to the direct impingement regions. 

V. Concluding Remarks 

A numerical simulation for "running wet" aircraft 
anti-icing systems was developed. The model includes 
breakup of the water film, which exists in regions of direct 
impingement, into individual rivulets. The wetness factor 
distribution resulting from the film breakup and the rivulet 
configuration on the surface were predicted in the 
numerical solution. The solid wall was modeled as a 
multi-layer structure and the anti-icing system used was of 
the thermal type utilizing hot air and/or electrical heating 
elements embedded within the layers. The mathematical 
formulation of the heat transfer process as well as details 
of the numerical solution procedure were presented. 

Experimental tests were recently conducted in the 
NASA Lewis Icing Research Tunnel to validate the current 
model. A detailed comparison with the numerical results 
was not possible at the time this manuscript was written 
since the data acquired were not reduced. However, similar 
trends were observed between the computer code 
predictions and the experimental results. Further detailed 
comparisons will be carried in the near future. 
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Table 1: Composite Wall Physical and 
Thermal Properties. 


Layer 

number 

Description 

Material 

Thickness 

(mm) 

Therm aJ 
Conductivity 
(W/m.K) 

1 

Substrate 

BSI 

2.20 

220 

2 

Inner 

Insulation 

Epoxy/ 

Glass 

1.30 

1.25 


Heater* 

Copper 

0.20 

340 

| 

Outer 

Insulation 

Epoxy/ 

Glass 

0.25 

1.25 

1 5 

Abrasion 

Shield 

Stainless 

Steel 

0.30 

50 


Heater turned off 
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Fig. 1: Rivulet and equivalent rectangular film models. 



Fig. 2: Type of grid points in the wall. 



x / chord 


Fig. 3: NACA 0012 airfoil geometry. 
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Wrap distance (m) Wrap distance (m) 

Fig. 10: Surface heat flux. Fig. 13: Surface water mass flow rate. 



Wrap distance (m) 


Fig. 11: Runback water film thickness. 



Fig. 12: Wetness factor. 



